subroutine Force(X, V, TM, F)
    implicit none

    REAL*8 X(3), V(3), F(3), TM  !in

    INTEGER body, coord
    REAL*8 r,rr,rrr,d(3)

    real*8 GM(11)
    COMMON /GM/GM

    REAL*8 RegBodyX(3,11), RegBodyV(3,11)
    INTEGER CenterBody
    COMMON /CenterBody/CenterBody

    logical(1) bUseRelative
    DATA bUseRelative/.FALSE./
!--------------------------------

    do coord = 1, 3
        F(coord) = 0.0D0
    end do

    call GetCoordVelRegBody(TM, RegBodyX, RegBodyV, CenterBody)

    if( CenterBody .EQ. 0) then ! bary center
        do body = 1, 11
            rr = 0.0d0
            do coord = 1, 3
                d(coord) = RegBodyX(coord, body) - X(coord)
                rr = rr + d(coord) * d(coord)
            end do
            r = dsqrt(rr)
            rrr = r*rr
            do coord = 1, 3
                F(coord) = F(coord) + GM(body)*d(coord)/rrr
            end do
        end do

        if(bUseRelative) then
            call RelativeForce(TM, X, V, GM, RegBodyX, RegBodyV, F)
        end if
    else
        rr = 0.0d0
        do coord = 1, 3
            rr = rr + X(coord) * X(coord)
        end do
        r = dsqrt(rr)
        rrr = r*rr
        do coord = 1, 3
            F(coord) = - GM(CenterBody) * X(coord) / rrr
        end do

        do body = 1, 11
            if( body .EQ. CenterBody ) cycle

            rr = 0.0d0
            do coord = 1, 3
                d(coord) = RegBodyX(coord, body) - X(coord)
                rr = rr + d(coord) * d(coord)
                !write (*,*) TM, body, coord, X(coord), RegBodyX(coord, body), d(coord), rr
            !    write (*,*) TM, body, coord
                !write (*, '(d29.20,1x,d29.20)') d(coord), rr
            end do
            r = dsqrt(rr)
            rrr = r*rr
            do coord = 1, 3
                F(coord) = F(coord) + GM(body) * d(coord) / rrr
                !write (*,*) TM, body, coord, X(coord), RegBodyX(coord, body), d(coord), rr,  F(coord)
            end do

            rr = 0.0d0
            do coord = 1, 3
                d(coord) = RegBodyX(coord, body)
                rr = rr + d(coord) * d(coord)
            end do
            r = dsqrt(rr)
            rrr = r*rr
            do coord = 1, 3
                F(coord) = F(coord) - GM(body) * d(coord) / rrr
            end do
        end do
    end if
end

! =========================================================

subroutine RelativeForce(TM, X, V, GM, RegBodyX, RegBodyV, Fres)
    implicit none
    REAL*8 TM, X(3), Fres(3), V(3)
    real*8 GM(11)

    INTEGER body, body2, coord
    REAL*8 r,rr,rrr,d(3)
    REAL*8 ForceAdd(3)

    REAL*8 RegBodyX(3,11), RegBodyV(3,11)
    REAL*8 Beta, Gamma, c, cc

    logical TESTOUT
    DATA     Beta/1.0e0/
    DATA    Gamma/1.0e0/

    DATA    TESTOUT/.FALSE./
    DATA    c /63239.7787e0/ ! ae/year
    !DATA    c /299792458.0e0/

    REAL*8 K1, K2, K3, K4, K5, K6, K7
    REAL*8 tmp

    real*8 kmAE
    common /kmAE/kmAE


!--------------------------------

    cc = c * c

    K3 = 0.0
    do coord = 1, 3
        K3 = K3 + V(coord) * V(coord) / cc
    end do
    K3 = K3 * Gamma

    do body = 1, 11
        rr = 0.0
        do coord = 1, 3
            d(coord) = RegBodyX(coord, body) - X(coord)
            rr = rr + d(coord) * d(coord)
            ForceAdd(coord) = 0.0
        end do

        K1 = 0.0
        K2 = 0.0
        K4 = 0.0
        K5 = 0.0
        K6 = 0.0
        K7 = 0.0
        do coord = 1, 3
            K4 = K4 + RegBodyV(coord, body) * RegBodyV(coord, body) / cc
            K5 = K5 + V(coord) * RegBodyV(coord, body)
            K6 = K6 - d(coord) * RegBodyX(coord, body)
        end do

        r = dsqrt(rr)
        rrr = r * rr

!        do coord = 1, 3
!        end do

        K6 = K6  / rr;
        K6 = K6 * K6;

        do body2 = 1, 11
            K1 = K1 + GM(body2) / r
            if( body .NE. body2) then
                K2 = K2 + GM(body2) / r
            end if
        end do

        K1 = (2.0) * (Beta + Gamma)/cc * K1
        K2 = (2.0 * Beta + 1.0)/cc * K2
        K4 = (1.0 + Gamma) * K4
        K5 = 2.0 * (1 + Gamma) / cc * K5;
        K6 = 1.5 * K6 / cc
        K7 = 0.5 / cc * K7

        do coord = 1, 3
            ForceAdd(coord) = ForceAdd(coord) + (- K1 - K2 + K3 + K4 - K5 - K6 + K7) * GM(body) * d(coord) / rrr
        end do
        if( TESTOUT .AND. body .EQ. 11) then
            !write (*,*) TM
            !write (*,'(d29.20)') K1, K2, K3, K4, K5, K6, K7
        end if
    end do

    do body = 1, 11
        tmp = 0.0
        do coord = 1, 3
            tmp = tmp - d(coord) * ((2.0 + 2.0 * Gamma) * X(coord) - (1.0 + 2.0 * Gamma) * RegBodyV(coord, body) )
        end do
        tmp = tmp / cc
        do coord = 1, 3
            ForceAdd(coord) = ForceAdd(coord) - tmp * GM(body) * d(coord) / rrr
        end do

        if( TESTOUT .AND. body .EQ. 11) then
            !write (*,'(d29.20)') tmp
        end if
    end do

    if (TESTOUT) then
        do coord = 1, 3
            Fres(coord) = Fres(coord) + ForceAdd(coord)
        !write(*,*) Fres(coord), ForceAdd(coord) , ForceAdd(coord)  / Fres(coord)
        end do
    !write (*,*) ''
    endif
end
